Wave breaking and particle jets in intense inhomogeneous charged beams 
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This work analyzes the dynamics of inhomogeneous, magnetically focused high-intensity beams 
of charged particles. While for homogeneous beams the whole system oscillates with a single fre- 
quency, any inhomogeneity leads to propagating transverse density waves which eventually result in 
a singular density build up, causing wave breaking and jet formation. The theory presented in this 
paper allows to analytically calculate the time at which the wave breaking takes place. It also gives 
a good estimate of the time necessary for the beam to relax into the final stationary state consisting 
of a cold core surrounded by a halo of highly energetic particles. 



It is well known that magnetically focused beams of 
charged particles can relax from non-stationary into sta- 
tionary flows with the associated particle evaporation [l[ . 
This is the case for homogeneous beams with initially 
mismatched envelopes flowing along the magnetic sym- 
metry axis of the focusing system. Gluckstern Q showed 
that initial oscillations of mismatched beams induce for- 
mation of large scale resonant islands [!, [1] beyond the 
beam border: beam particles are captured by the res- 
onant islands resulting in emittance growth and relax- 
ation. A closely related question concerns the mecha- 
nism of beam relaxation and the associated emittance 
growth when the beam is not homogeneous. On general 
grounds of energy conservation one again concludes that 
beam relaxation takes place as the coherent fluctuations 
of beam inhomogeneities are converted into microscopic 
kinetic and field energy [JlH, d,0|- However, unlike in the 
former case for which the specific resonant mechanism is 
well understood, for inhomogeneous systems a more de- 
tailed description of the processes involved must still be 
explored. This is the goal of the present letter. 

The interest surrounding a better understanding of the 
dynamics of inhomogeneous beams is due to the fact 
that one can hardly design experimental devices capa- 
ble of generating fully matched beams at the entrance of 
a transport system [8|. While the azimuthal symmetry 
with respect to the beam axis is feasible in the case of 
focusing solenoidal magnetic fields, envelope matching — 
to avoid the radial oscillations — is significantly harder 
to achieve, while a perfect homogeneity is practically im- 
possible. 

Given all these facts, the purpose of the present work 
is to investigate the mechanisms leading to the decay of 
density inhomogeneities as the system relaxes into its fi- 
nal stationary state. We find that the relaxation comes 
about as a consequence of breaking of density waves fol- 
lowed by ejection of fast particle jets. Jets are formed by 
particles moving in-phase with the macroscopic density 
fluctuations. They draw their energy from the propagat- 



ing wave fronts and convert it into microscopic kinetic 
energy. This process is very similar to the breaking of 
gravitational surface waves. The jet can then be com- 
pared to a broken crest of a gravitational wave surfing 
down the wave front. We stress that the wave breaking 
mechanism analyzed in this letter is very different from 
the Gluckstern resonances which were found to be the 
driving force behind the omittance growth in transversely 
oscillating homogeneous particle beams. For strongly in- 
homogeneous beams, we find that it is the wave breaking 
and jet production which are the primary mechanisms 
responsible for the beam relaxation. 

We consider solenoidal focusing of space-charge dom- 
inated beams propagating along the transport axis, de- 
fined as the z axis, of our reference frame. The beam 
is initially cold with vanishing emittance, and is az- 
imuthally symmetric around the z axis. Since the number 
of constituent particles is very large, the beam dynamics 
is governed by the azimuthally equation and collective 
effects are dominant. Prior to the appearance of density 
singularities, the original Vlasov formalism can be sim- 
plified to a cold fluid description for which Lagrangian 
coordinates are particularly appropriate. In these coordi- 
nates, the transverse radial position r of a beam element 
is governed by 
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The prime indicates derivative with respect to the longi- 
tudinal z coordinate which for convenience we shall also 
refer to as "time" , and angular momentum in the Larmor 
frame is taken to be zero for each particle. The focus- 
ing factor is k = (qB /2~fmf3c 2 ) 2 , where B is the axial, 
constant, focusing magnetic field; Q(r ) = KN(r )/N t , 
is the measure of the charge contained between the ori- 
gin at r = and the initial position r(z = 0) = ro, 
N t is the total number of beam particles per unit axial 
length, N(ro) is the number of particles up to ro, and 
K = N t q 2 /"f 3 mf3 2 c 2 is the beam perveance. q and m 
denote the beam particle charge and mass, respectively; 
7 = (1 — /3 2 ) -1 / 2 is the relativistic factor where /3 = v z /c 
and v z is the constant axial beam velocity and c is the 
speed of light. Note that r n i s in fact the Lagrange coordi- 
nate of the fluid element [l3| which means that as long as 
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the fluid description remains valid, the amount of charge 
seen by the fluid element inside the region < r < r(z) 
remains unaltered at Q(ro), independent of time z. This 
is of fundamental importance since from the Gauss law 
this is the charge that exerts the force on the fluid ele- 
ment. In this letter we will consider the beams starting 
from a static initial condition, r'(0) = 0. The formal so- 
lution to the fluid dynamics Eq. ([T]) is r = r(z,ro). This 
can be calculated explicitly using the Lindstedt-Poincare 
perturbation theory. For small amplitude fluctuations 
around the stable equilibrium r eq (r ) = y/ Q(ro)/ n we 
obtain 

r(r ,z) = r eq {l + A/r eq cos(uz) + (1/3) {A/r eq f x 

[2 + cos(wz)] sin 2 (wz/2) + 0[{A/r eq f]} , (2) 

where A(ro) = ro — r eq is the amplitude of oscillations, 
w(ro) = ojq + y/~KA 2 /(6\/2 r^ q ) is the renormalized tq- 
dependent frequency, and uiq = \/~2k is the unperturbed 
frequency. We stress that as long as the fluid picture 
applies, all the information about the temporal evolu- 
tion of the beam is contained in Eq. @. For example, 
the time evolution of the beam density can be obtained 
as follows. For a beam of initial cross sectional density 
no(r), the amount of charge SQ between two concen- 
tric circles of radii ro and tq + Sro (<5ro small) is SQ = 
2TTr 5ro(K/N t )no(ro)- Since this charge is conserved, 
SQ = 2nrSr (K/N t )n{r) = 2nr{dr /dr a ){K/N t ) n{r)Sr , 
the transverse beam density at any future time z is, there- 
fore, 

n(r) = n (r ) (r Q /r) {dr/dr^ 1 \ ro=ro{r zy (3) 

For a given position r and an axial coordinate z, the 
initial position ro of a beam element can be uniquely de- 
termined as the inverse function tq — fo(r, z) of Eq. 
This ceases to be the case if dr/dro — > and r (r, z) 
becomes multivalued. If this happens, the density will 
diverge and the fluid picture will break down. All these 
features, if present, would be indicative of a wave break- 
ing phenomenon. Needless to say that presence of wave 
breaking in charged particle beams would be of consider- 
able interest and practical importance. Breaking might 
be responsible for conversion of energy from macroscopic 
fluid modes into microscopic kinetic activity. 

We start our analysis by considering the compressibil- 
ity factor dr/drQ, which can be obtained exactly by nu- 
merically integrating two nearby trajectories of Eq. ([I]) or 
approximately by differentiating Eq. @. To be specific 
we write the initial cross sectional beam density at z = 
in a general parabolic form no(ro) = pu + XPi( r o), where 
the inhomogeneity parameter < x < 1) Ph = N t /(irr%), 
and Pi(ro) = ph(2r\/r\ — 1). r\> is the beam radius. 
Note that the integral of the inhomogeneous contribution 
Pi(r ) is zero. To suppress the effects arising from the 
pure envelope oscillations — Gluckstern resonances - 
we fix n = y K/k, so that the beam radius is unaltered 
for as long as the fluid picture, Eq. (fTJ, remains valid. 
We have also performed calculations for rms matched 



beams and find that they behave qualitatively the same 
way. Fig. [1] shows the typical time evolution of the com- 
pressibility dr/drg obtained both numerically using Eq. 
(fl]) and analytically using $2$. Since the amplitudes of 
oscillations of all fluid elements about the points of their 
equilibria are small, \A/r eq \ <C 1, even for large values of 
X, the agreement between the numerical and the pertur- 
bative solution is found to be very good, see Fig. [TJ The 
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FIG. 1: Time evolution of the compressibility factor dr/dro 
for x — 0.6 and ro = OAvb- dr/dro = indicates singular 
density. Symbols represent the compressibility obtained from 
the numerical solution of Eq. |T} , whereas the solid line is de- 
rived using the perturbative solution, Eq. {5). z is measured 
in units of k _1//2 . Position ro = 0.4 is chosen because this 
is the approximate Lagrange coordinate at which density is 
first found to diverge for \ — 0-6. 

compressibility factor exhibits a fast oscillatory motion 
accompanied by a slow secular growth. This means that 
given enough time, the compressibility will always be- 
come zero for any finite value of x, resulting in a density 
divergence. 

To further explore the significance of the diverging den- 
sity, we have performed fully self-consistent N-particle 
simulations. For a system in which particles interact 
by an infinite range unscreened Coulomb potential, the 
time of collision diverges and the mean-field Vlasov de- 
scription becomes exact [14] , Thus, in order to simulate 
the Nt — ► oo limit — in which the thermalization time 
due to binary collisions is infinite — each particle can 
be taken to interact only with the mean-field produced 
by all the other particles. Taking advantage of the az- 
imuthal beam symmetry and the Gauss law, a particle 
located at a position r experiences a field generated by 
the particles within a circle of radius r 0, [lj| . We stress 
that in the Nt — > oo limit this is exact at any finite 
time scale. Within the simulation, the trajectory of each 
particle is, therefore, also governed by Eq. |T]) — un- 
like the fluid elements, however, particles are allowed to 
bypass one another. This method avoids the thermal- 
ization effects associated with the binary collisions and 
significantly speeds up the simulations. It allows us to 
accurately simulate the collective effects dominant for all 
time scales when Nt — > oo even with a relatively small 
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number of particles, N t ~ 10000. 

In Fig. [2]we display the particle phase-space (r, v = r') 
for x — 0-6. The first panel (a) shows the initial distribu- 
tion at z = — all particles are still. In panel (b), after 
various propagating wave cycles, the system is about to 
build up an infinite density; velocity is still a single val- 
ued function of the space coordinate but a singularity 
(cusp) is forming. The third panel (c) shows the sys- 
tem at a time slightly larger than the first wave break- 
ing time. Velocity ceases to be a single valued function 
of the space coordinate — while some of the particles 
go through the wave others do not. This latter class 
of particles is accelerated by the wave front and forms 
a thin azimuthally symmetric jet or finger seen in the 
figure. High energy jet particles can reach far outside 
the beam core and may be very detrimental to the beam 
transport. The process shown in panel (c) repeats itself 
many times, see panel (d), as the system evolves toward 
a final stationary state and the previously unoccupied 
extensions of the phase-space are gradually filled with 
particles whose velocities are considerably larger than ve- 
locities in the beam core. After some time a stationary 
state is reached in which the beam separates into a cold 
dense core and a hot and extended halo of ejected par- 
ticles, panel (e). Time evolution of the emittance 
e = 2 V < r 2 >< v 2 > — < rv > 2 , where < > denotes an 
average over particles, is shown in the last panel (f). 
At the wave breaking emittance suffers a sharp rise, fol- 
lowed by a rapid relaxation to the final stationary state 
in which large amplitude fluctuations subside. We note 
that beam relaxation is closely connected to phase-space 
filamentation. Phase-space filamentation takes place af- 
ter particles are ejected from the beam core by surfing 
on the charge density waves. Once outside the core, 
particles experience the time dependent nonlinear forces 
and undergo all the complicated mixing dynamics with 
subsequent filamentation. This leads to final irreversible 
emittance growth. The directed emittance growth seen 
prior to relaxation is reversible. Since the beam radius 
rf, is matched — envelope oscillations are small — the 
contribution of the Gluckstern resonant mechanism to 
the emittance growth and beam relaxation is, at most, 
marginal. The dominant mechanism is the singular build 
up of density followed by the wave breaking and jet pro- 
duction. The time of the first wave breaking depicted in 
the panel (c) of Fig. [5] agrees well with the time when 
the compressibility factor dr/dr^ obtained from the La- 
grange fluid equations goes to zero, Fig. [TJ Our next 
goal is then to precisely calculate the instant at which 
the wave breaking takes place. 

We first note that for an inhomogeneous density pro- 
file, each fluid element oscillates with a different fre- 
quency — rigid oscillations are possible only when the 
density profile across the beam cross section is homoge- 
neous. Thus, nearby fluid elements will oscillate around 
their points of equilibria, slowly moving out-of-phase. 
This motion results in transverse density waves propa- 
gating across the beam. At some point, however, two 
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FIG. 2: Time evolution in the phase-space (r, v) for \ = 0.6. 
r is measured in units of (K/k) 1 ^ 2 , v in units of K 1 ^ 2 , and z in 
units of jt -1//2 : panel (a) - initial condition at z = 0; panel (b) 
- after many wave cycles, but just before wave breaking; panel 
(c) - just after wave breaking, an azimuthal jet is expanding 
over the phase-space; panel (d) - while the first jet moves 
in the phase space new jets are being ejected; panel (e) - a 
phase portrait of a relaxed state; panel (f ) - emittance growth 
leading to a final relaxation. Note that emittance saturates 
soon after the wave breaking (wb). 



nearby fluid elements will overlap one another leading to 
a singular build up of density. When this happens the 
fluid picture will lose its validity and will have to be re- 
placed by the full kinetic description given in terms of 
the Vlasov equation. The wave breaking occurs when 
the separation between any two fluid elements vanishes, 
r(ro + <5ro, z) — r(Yo, z) — > 0, for some value of r$. This 
is precisely equivalent to our condition for the appear- 
ance of a singular density, dr/dro — > 0. Considering 
only the term linear in amplitude of Eq. ([2]) , we see that 
Sr = Sr eq + cos(ujz)SA — Az sm(ojz)5cu. Neglecting the 
purely oscillatory term, as compared to the secular one, 
the time of breaking is found to be 
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As expected, the breaking will always occur whenever 
duj/dr 7^ 0. This is the case for all inhomogeneous par- 
ticle beams. Unlike other systems in which one must 
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have strong enough electric fields [16|], here any sort of 
inhomogeneity leads to the wave breaking — one just has 
to wait sufficiently long. As soon as the wave breaking 
takes place, particles with the same velocity as the den- 
sity wave will be captured by the wave and surf down its 
front gaining kinetic energy (Fig. [2]). Since at the wave 
breaking position ro/r b < 1, minimization in Eq. ^ can 
be performed perturbatively in this parameter, yielding 

_{3\ 1/2 a 3 (4yTf y + X~l) 
^ W (V3- a ) 2 (VT^ + X-l)' ^ 

where a = (1 + 2 \Jl - x ~ X) 1 ^ 2 - 

For small values of X the nrs t wave breaking will hap- 
pen after a long transitory period; z w b » 81%/2/x 3 as 
X — > 0. The subsequent breaks, however, occur on a 
much shorter time scale ~ l/u>o, as can be seen from Fig. 
[TJ Therefore, z wb should also give us a good estimate of 
the relaxation time for the entire dynamics. In Fig. [3] 
we compare the wave breaking time obtained using the 
N particle dynamics simulation described above with z w b 
given by Eq. ([5]). The figure reveals an amazingly good 
agreement between the two results. In the same plot we 
also show the relaxation time — defined as the time when 
the emittance first reaches its plateau value, see Fig. [5] 
(f). As expected, for smaller values of x the time of re- 
laxation follows closely the wave breaking time. This is 
because the phase mixing and the jet production occur 
on a much shorter time scale, ~ 1/wq, than z w b- For 
larger values of x the two time scales, however, become 
comparable. This results in a deviation between the two 
data sets — circles and crosses — observed in Fig. [3] 



To summarize, we have investigated the dynamics of 
space-charge dominated beams |10| with inhomogeneous 
density profiles. Using Lagrangian coordinates, we were 
able to derive a very accurate analytical expression Eq. 
([2]) which describes very well the dynamics of beam par- 
ticles, up to the wave breaking time. The fluid picture 
loses its validity when the propagating wave fronts result 
in a singular build up of density. At this point the crest 
of the propagating wave will break off producing an az- 
imuthally symmetric jet of particles accelerated by the 
wave front. This process will repeat itself many times 
leading to a final stationary state in which the beam sep- 
arates into a cold core surrounded by a halo of highly 
energetic particles. The theory presented in this paper 
allows to calculate precisely the time at which the wave 
breaking will take place. It also gives a very good es- 
timate for the time of relaxation to the final stationary 
state. Unlike other systems in which the wave break- 
ing occurs only when thresholds on driving fields are ex- 
ceeded [16|. inhomogeneous beams are found to be always 
unstable [l7[ and the wave breaking is unavoidable. 



Wave breaking is not the only mechanism which leads 
to the relaxation of initially non-stationary intense par- 
ticle beams. It is well known that oscillations of mis- 
matched envelopes can be damped by the Gluckstern 
resonances. However, in practice, inhomogeneities are 
much harder to suppress than the envelope mismatches 
8]. In these cases the wave breaking described in the 
present letter will be the dominant mechanism by which 
a system reaches its final stationary state. 




FIG. 3: Comparison of the predicted time of the first wave 
breaking Eq.([5} (solid line), with the result of dynamics sim- 
ulations (circles). Time is measured in units of k _1//2 . A very 
good agreement between the theory and the simulations ex- 
tends all the way to x = 0.8. Crosses show the relaxation 
time — when emittance saturates. Note that up to x = 0.6 
the wave breaking time and the relaxation times are almost 
identical, see the discussion in the text. 



The theory presented here describes beams with van- 
ishing initial emittances. One example of this are the 
crystalline beams for which the initial emittance is sup- 
pressed by a series of dissipative cooling procedures [18| • 
We expect, however, that the theory will remain valid 
also for beams of initial finite emittances as long as the 
thermal lengt h is sma ll, v th z wb = e th z wb /2r b < A ~ x r b, 
where Vth — V< v 2 > is the characteristic thermal veloc- 
ity, A is the characteristic amplitude of particle oscilla- 
tions, and the value of the wave breaking time is given in 
Fig. ([3]). Since in general z wb ^> 1/wo and x < 1) the con- 
dition requires that the thermal velocity be smaller than 
the macroscopic velocity uqA and that the beam be 
space charge dominated e <C KJ y/n. 
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